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Abstract 

We obtain an analytical solution for the pressure-transient behavior of a vertically hy¬ 
draulic fractured well in a heterogeneous reservoir. The heterogeneity of the reservoir is mod¬ 
eled by using the concept of fractal geometry. Such reservoirs are called fractal reservoirs. 
According to the theory of fractional calculus, a temporal fractional derivative is applied to 
incorporate the memory properties of the fractal reservoir. The effect of different parameters 
on the computed wellbore pressure is fully investigated by various synthetic examples. 

Keywords: vertically hydraulic fractured well; fractal geometry; fractal reservoir; fractional 
derivatives. 


1 Introduction 


Hydraulic fracturing plays an important role in improving the productivity of damaged wells and 
wells producing. The vertical plane fracture is created by injecting fluid into the formation and 
then filling with proponing agents, such as propants, to prevent closure. In practical terms, two 
types of fractured well are considered: infinite (high) or finite (low) conductivity vertical fracture. 
In case of infinite conductivity fracture, it is assumed that the fluid flows along the fracture without 
any pressure drop. Finite conductivity fracture occurs when the pressure drop along the fracture 
plane is not negligible. 

The classical diffusion equation has been used to explain the pressure responses of a well in a 
reservoir, which is assumed to be homogenous at all scales. However, recent studies show that the 
homogeneity assumption is not valid in most cases [ljj7|. Due to this fact, fractal geometry has 
been used as an effective tool to describe the heterogeneities of these reservoirs, which are called 
fractal reservoirs [&-11 . Since the diffusion process of fractal reservoirs is history dependent, and 


the anomalous diffusion properties of fractal reservoirs cannot be fully described by the fractal 
model, the concept of fractional derivative has been used to incorporate the memory of the fluid 
flow 12-15 . Our main aim here is to analyze the pressure behavior of a well with an infinite 
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conductivity vertical fracture in a fractal reservoir. An infinite radial system is considered in order 
to analyze the effects of different parameters on the well response. 

The paper is organized as follows. After this brief introduction, the mathematical model 
is formulated in Section [2] A summary of the nomenclature used appears in Appendix [A] The 
analytical solution to the model is provided in Section[3] Section [4] discusses how the well responses 
to different reservoir parameters. The main conclusions of our study are given in Section [5] 


2 Model description 

A schematic diagram of a vertically hydraulic fractured well is shown in Figure [l] Figure [2] shows 
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Figure 1: Geometry of a vertically hydraulic fractured well, 
the geometry of flow lines near the fractured well. Before discussing the mathematical model of 




Figure 2: Linear and pseudo-radial flow regimes near an infinite conductivity fracture. 


transport process, we define three variables: 


• the dimensionless pressure 


Pd 
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qp 


• the dimensionless time 
































• the dimensionless radius 


See Appendix [A] for the description of all the quantities involved. The best known and most 
useful model to describe the pressure behavior of fractal reservoirs was firstly proposed by Metzlcr 
et al. [l]. Similarly to Camacho-Velazquez et al. [4], from here on we call generalized diffusion 
equation to the fractal-fractional diffusion (FFD) equation 


1 d 2 p D /3 dp D _ d 1 pp 

r D dr D r ^ +1 dr D dt D 7 ’ 


where /3 = d m f — 0—1. The parameter d m f denotes the mass fractal dimension, while 0 represents 
the conductivity index. Mass fractal dimension is responsible for the reservoir structure, and the 
conductivity index explains the diffusion process in the reservoir. The Caputo fractional order 
derivative is used to introduce d^pp / dtp 1 '■ 
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where 7 £ M + , [ 7 ] = m £ Z + , and T represents the Gamma function, that is, 
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The order 7 of the fractional derivative is related to the conductivity index by 7 = 2/(2 + 0). 


3 Analytical solution 


It is assumed that the pressure distribution of the reservoir is uniform and constant at initial time: 


Vd{t l >,0) = 0. (4) 

To obtain the line source solution of equation ([!]), we take r w —> 0 + . The inner boundary condition 
without wellbore storage and skin effects can then be written as 


rD 


0 dp D {r Dl t D ) 
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= - 1 . 


r D ->-0+ 


( 5 ) 


The outer boundary condition for the infinite reservoir is given by 


lim PD{r Dl t D ) = 0 . 

rr>—>00 


Taking the Laplace transform to both sides of ([!]), and then using Q, we obtain that 

1 d 2 p D /3 dp D 7 _ 
r e D dr% + r e D +1 dr D * Pd ’ 
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where s is the Laplace transform variable. The dependent variable pp denotes the Laplace trans¬ 
form of pd, and is a function of rp and s. In the Laplace space, the inner boundary condition 
takes the form 


p dp D (r D ,s ) 
dr D 


r*£)—>-0+ 
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( 8 ) 


while the outer boundary condition is given by 


lim pp{r Dl s) = 0. (9) 

fD-yoo 


3 













Using the substitutions pn = r d^ 1 P^ 2 W and x = (2s 7 / 2 rc/ 2+e )/ 2 )/(2 + 0), and after a slight 
manipulation, we conclude that © is equivalent to 


,d 2 W dW 


r „ , + x—-(x + v~)W = 0, 


dx 2 dx 

where v = ^f§. Equation (10) is Bessel’s equation, which has the general solution 

W(x, s ) = AI„(x) + BK v {x). 


Thus, the dimensionless pressure function in Laplace space can be written as 


Pd{t D ,s) = t\d (1 13)12 
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Application of the outer boundary condition Q to (12) yields A = 0. Therefore, (12) reduces to 

Pd(td,s) = Br D ^-~ p)/2 K v ^^- 0 rz} (2+e)/2 ^ • (13) 

Based on (13), the inner boundary condition (| 8 j) can be written as 

— B lim s^ 2 r D )1+ ^ +e) / 2 K v _ l ( 2 ^- rD )2+e )/ 2 \=--. 

r D -)-0+ \2 + 0 / S 


(14) 


Equation (14) can be simplified by using the formula 


K v (x) 


I» /2V 


v > 0 , 


(15) 
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valid for small arguments 0 < x <C \j\ + v. Indeed, having in mind that AT„(x) = K_ v (x), and 
by making use of (15), we reduce (14) to the following expression: 

l-v 
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Equation (13) can be written as 


T ( l -^)(2 + ( 9 ) s 1 ^ 7 / 2 
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Pd(xd,Ud,s) = B (x D - aif + (y D - a 2 ) 2 ^ 

/ p q 7/ 2 / /-\ (2+ e )/2\ 

[j^iy/ixD-a^ + iyn-a^j j. (18) 


x K„ 


We obtain the pressure drop function, in a fractal reservoir with a hydraulic fracture across the 
well, by integrating (13) with respect to a\ from —1/2 to 1/2: 

(l-/3)/2 
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Since the computation of the pressure along the fracture is favorable, it can be assumed that 

(2+0)/2N 


Ud = ol 2 - Thus, the pressure drop function (|19| reduces to 
ri /2 
1 - 1/2 


rl/2 / , -\ (1-^/2 (2 S 7/ 2 / /-NCi-KM/z \ 

Pd{x d ,0iS) = B J ^ y\J (x D - Oiif j K v ( \j (x D - a^ 2 J J da i, (20) 

where |:ed| < 1/2. By changing variables, and after further manipulations, (20) can be written as 
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with a = (4 — + 0)/(2 + 9). The integrals in (21) can be computed by the following formula: 
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where v ^ 'L and 1 -F 2 represents the generalized hypergeometric function. For the Euclidean 
model, that is, the particular case v = 0 and a = 1, the presented integrals can be evaluated by 


J K 0 (z)dz = y {K 0 (z)L_ 1 (z ) + A'i(2)i 0 (^)) 


or 


J K 0 (z)dz = ^ (AT 0 (z)Ai( 2 ) + Ki(z)L 0 {z)) + zK 0 (z), 
where L_i( 2 ), Lq{z) and Li(z) are modified Struve functions given by 
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The inverse transform operation is easily carried out by using the Gaver-Stehfest algorithm 16 


4 Model responses 


To show the pressure-transient behavior of a fractured well, we use (21) to compute the wellbore 


pressure response at Xd = 0 (uniform flux solution), without wellbore storage and skin effects. 
The well intercepts the center of a single vertical fracture plane. According to the uniform flux 
fracture assumption, the flow rate per unit of fracture surface is constant along the fracture length. 
Figure [3] shows the different behavior of a well response in a homogenous and fractal reservoir. 
Generally, the fracture structure works as a sink, which enforces the fluid to go towards the 
fracture with hyperbolical flow geometry. The geometry of fracture enforces the flow lines to be 
perpendicular to the fracture plane and the pressure-transient response defines a linear flow in the 
reservoir. After the linear flow regime, when the effect of fluid flow reaches the two ends of the 
fracture, the geometry of flow can be identified, for a small time, by a hyperbolical flow regime. 
Ultimately, the characteristic radial flow regime is observed from the well response. The linear flow 
regime can be identified by two straight-lines of pressure and its logarithmic derivative at early 
times (see Figure [3]). As can be seen from Figure [3j these two parallel straight-lines have a half 
slope during the linear flow regime in a conventional reservoir. The specialized analysis (Cartesian 
coordinates) on Figure [4] shows that with a plot of the pressure change versus the square root 
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of the elapsed time, the linear flow can be identified with a straight-line intercepting the origin. 
The results of Figure [4] show that the slope of the straight-line is decreasing with the increasing 
of the parameter v. In other words, in a more complex reservoir, described by a larger value of 
the conductivity index 9 and a smaller value of the mass fractal dimension d m f , the linear flow is 
identified with a smaller slope of the straight-line in the specialized plot. Figure [3] indicates that 
the radial flow regime (at late times) is identified by the straight-line shape of the logarithmic 
derivative in log-log scale and that the reservoirs with more complex structures (larger values of 
v) have the straight-lines with larger slopes. Moreover, Figure [3] shows that in a fractal reservoir, 
at late times during the radial flow regime, the pressure response and its logarithmic derivative 
are parallel straight-lines. This fact indicates that the pressure and its logarithmic derivative can 
be expressed by two power-law functions with the same power. 



Figure 3: Responses for a fractured well without wellbore storage and skin effects. Log-log analysis 
of pd and dpn/dln(tp) versus to with d m f = {2, 1.95, 1.9, 1.85, 1.8}, 9 = {0, 0.1, 0.3, 0.5, 0.7} 
and 7 = {1, 0.9524, 0.8696, 0.8000, 0.7407}. 

The reservoir rock porosity is a statistical property of the system, and can be expressed as 
</>(r) = <fi r (r/L r ) dmf 2 , where <p r is a reference porosity value at the reference length L r . Partic¬ 
ularly, the reference length is the wellbore radius, and the porosity of a fractal reservoir can be 
written as </>(r) = <t> w {'r/r w ) drri:f 2 . The reference porosity value can be assumed to be equal to the 
porosities obtained from log and core data. However, based on log and core porosities, the reason¬ 
able oil-in-place estimates indicate that the porosity variations are not big, and can be considered 
equal to the constant reference porosity. It can be therefore immediately concluded that the mass 
fractal dimension d m f is equal to the Euclidean dimension, that is, d m f = 2. On the other hand, 
because k{r) = k r (r/L r ) drnf 9 2 and, particularly, k{r) = k w (r / r w ) dmf 9 2 , the deficiencies of 
average permeability calculations indicate that the dynamical properties of the system have an 
important role in the flowing fluid through the reservoir and in all stages of production. Since 
d m f = 2, the permeability variations can be expressed by k(r) = k w (r/r w )~ 9 . So, it requires 
investigating the effect of the conductivity index 9 variations on the well response. As can be seen 
from Figure [5] by increasing the 9 values, the pressure responses and their logarithmic derivatives 
increase. For a constant-rate production, more pressure drops in reservoirs with a more compli¬ 
cated diffusion process (larger 9) indicates that the diffusion is slower in these types of reservoirs. 
Therefore, application of the conventional model solutions to such reservoirs may lead the analyst 
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Figure 4: Specialized analysis of po versus y/to with d m f = {2, 1.95, 1.9, 1.85, 1.8}, 6 = 
{0, 0.1, 0.3, 0.5, 0.7} and 7 = {1, 0.9524, 0.8696, 0.8000, 0.7407}. 


to a wrong interpretation. 


10 


-O 10 

o 

'E 


’ 10 “ 


10 


m 10"' 


E 

Q 10 



- Pressure drop, nu=0 

-Logarithmic derivative, nu=0 

- Pressure drop, nu=0.0476 

- Logarithmic derivative, nu=0.0476 

-Pressure drop, nu=0.1304 

-Logarithmic derivative, nu=0.1304 

- Pressure drop, nu=0.2000 

-Logarithmic derivative, nu=0.2000 

-Pressure drop, nu=0.2593 

-Logarithmic derivative, nu=0.2593 


10 


10 


10 


10 10 10 
Dimensionless time 


10 


10 


10 


Figure 5: Responses for a fractured well without wellbore storage and skin effects. Log- 
log analysis of pn and dpo/dln(to) versus to with d m f = 2, 9 = {0, 0.1, 0.3, 0.5, 0.7} and 
7 = {1, 0.9524, 0.8696, 0.8000, 0.7407}. 
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5 Concluding remarks 

We used fractional calculus to model an infinite conductivity vertically fractured well in a fractal 
reservoir with more realistic behavior than the conventional (homogeneous) reservoirs. An ana¬ 
lytical solution was obtained, allowing us to analyze and interpret the well test data taken from 
these types of wells and reservoirs. On the basis of our results, the simulated examples and the 
discussion presented, the following conclusions can be expressed: 

• the Euclidean interpretation model cannot be used to analyze the pressure-transient behavior 
of fractal reservoirs; 

• the obtained solution can be used as a practical tool to analyze the effects of dynamical 
properties of the system; 

• according to the shortcomings of the average permeability estimates, the analytical solution 
can be used to determine the conductivity index parameter 6 and then the permeability 
variations. 


A Nomenclature 


Symbol 

Meaning 

Units 

c 

compressibility 

vol/vol/atm 

dmf 

mass fractal dimension 


1 F 2 

generalized lrypergeometric function 


h 

reservoir thickness 

cm 

h 

modified Bessel function of the first kind of order v 


k 

permeability 

darcy 

k w 

permeability at the wellbore 

darcy 

K v 

modified Bessel function of the second kind of order v 


L f 

fracture half-length 

cm 

K 

modified Struve function of order v 


V 

pressure 

atm 

PD 

dimensionless pressure 


Pi 

initial pressure 

atm 

1 

production rate 

cc/sec 

r 

radial distance 

cm 

I'D 

dimensionless radius 


r w 

wellbore radius 

cm 

T wD 

dimensionless wellbore radius 


t 

time 

sec 

tD 

dimensionless time 


S 

Laplace image space variable 


7 

fractional derivative order 


r 

Gamma function 


e 

conductivity index 


p 

viscosity 

cP 

0 

porosity (or void fraction) 


4*w 

porosity at the wellbore 
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